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1. Introduction 


The direct numerical simulation of turbulent flows at the high Reynolds numbers encoun- 
tered in problems of technological importance is all but impossible as a result of the wide 
range of scales that are present. Consequently, the solutions to such problems must invari- 
ably be based on some form of turbulence modeling. Traditional turbulence models based 
on Reynolds averages have had only limited success since the large scales of the turbulence 
(which contain most of the energy) are highly dependent on the geometry of the flow being 
considered. Experience has indicated that such models usually break down when a variety 
of turbulent flows are considered (Lumley 1983). The small scales are more universal in 
character, and serve mainly as a source for dissipation. Hence, it can be argued that a 
better understanding of turbulent flows could be achieved if just the small scales are mod- 
eled while the large scales are calculated (Deardorff 1970). This is the fundamental idea 
behind large-eddy simulations. 

During the past decade, considerable progress has been made in the large-eddy simu- 
lation of incompressible turbulent flows. This effort has shed new light on the physics of 
turbulence. The earliest work relied heavily on the use of the Reynolds averaging assump- 
tion to eliminate the Leonard and cross stresses while the Reynolds stresses were computed 
using the Smagorinsky model (Deardorff 1970, Leonard 1974, Reynolds 1976). More recent 
large-eddy simulations have been based on the direct calculation of the Leonard stresses 
with models provided for the cross and Reynolds subgrid-scale stresses in order to enhance 
the numerical accuracy (see Biringen and Reynolds 1981, Bardina Ferziger and Reynolds 
1983). However, among these models, only the particular Bardina, Ferziger and Reynolds 
(1983) model, with a Bardina constant of 1.0, satisfies the important physical constraint 
of Galilean invariance (Speziale 1985). The underlying physical concepts, fundamental nu- 
merical algorithms, and comprehensive historical data behind the recent field of large-eddy 
simulation have been presented in articles by Schumann (1975), Voke and Collins (1983) 
and Rogallo and Moin (1984). 

Despite the intensive research effort that has been devoted to the large-eddy simulation 
of incompressible flows as outlined above, it appears that no large-eddy simulation of a 
compressible turbulent flow has yet been attempted. Of course, such work could have 
important technological applications in the analysis of turbulent supersonic flows where 
shock waves are generated and in turbulent flows in combustion chambers. The prerequisite 
for carrying out such computations is the development of suitable subgrid-scale models for 
compressible turbulent flows. With the exception of the recent work of Yoshizawa (1986) 
few, if any, studies along these lines appear to have been published. The subgrid-scale 
models of Yoshizawa are only suitable for slightly compressible turbulent flows since they 
made use of an asymptotic expansion about an incompressible state. Hence, there is a 
need for subgrid-scale models that can be applied to the study of turbulent flows which 
are strongly compressible (i.e. for Mach numbers M > 1 where shock waves can occur). 
This forms the raison d’etre for the present study. 
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In this paper, subgrid-scale models will be developed for the closure of the Favre- 
filtered Navier-Stokes and energy equations. The compressible subgrid-scale stress model 
that will be obtained in Section 2 reduces to the linear combination model of Bardina, 
Ferziger and Reynolds (1983) in the incompressible limit. Likewise, the subgrid-scale heat 
flux model that will be obtained herein consists of an analogous linear combination of 
scale similarity and gradient transport terms. The dimensionless constants which appear 
in these subgrid-scale models will be arrived at by correlating with the results of direct 
numerical simulations of compressible isotropic turbulence. A more detailed comparison 
of computed and modeled subgrid-scale fields will be presented along with a discussion of 
the prospects for future research. 


2. Subgrid-Scale Models for Compressible Turbulence 


The compressible turbulent flow of an ideal gas will be considered. Such flows are governed 
by the continuity, momentum and energy equations given by (Batchelor 1967) 


dp d[pVk) 

dt dx k 
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respectively, where p is the mass density, v is the velocity vector, p is the thermodynamic 
pressure, p is the dynamic viscosity, h is the enthalpy, T is the absolute temperature, and 
k is the thermal conductivity. The viscous stress a k i and the viscous dissipation $ are 
respectively defined by 


<?ki 
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2 dv i c , ,dv k dvi 
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dxi dx k dxi 


( 4 ) 

( 5 ) 


The Einstein summation convention applies to repeated indices and the effects of external 
body forces have been neglected for simplicity. Equations (l)-(3) must be supplemented 
with the equations of state for an ideal gas which are as follows: 


P = pRT, h = C P T (6) 

where R is the ideal gas constant and C p is the specific heat at constant pressure. Likewise, 
the dependence of the viscosity and thermal conductivity on the temperature must be 
provided (i.e., relationships of the form p = p(T) and k = k(T) are needed and these 
depend on the gas under study). 
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( 7 ) 


Any flow variable 7 can be filtered in the following manner: 

T(x) = / G(x - z, A) 7(z)d 3 z 

where G is a filter function, A is the computational mesh size, and D is the domain of the 
fluid. The filter function G can either be an infinitely differentiable function of bounded 
support in an arbitrary domain, or a Gaussian distribution in a periodic domain It is 
normalized by requiring that 

f G(x— z, A)d 3 z = 1 . (8) 

Jd 

It follows that in the limit as the computational mesh size goes to zero, (7) becomes a 
Dirac delta sequence, i.e. 

lira G{x — z, A) 7 (z )d 3 z = J^ & (x — z) 7 (z)d 3 z = 7 (x) (9) 

where d(x — z) is the Dirac delta function (Arfken 1970) As a direct result of the Riemann- 
Lebesgue theorem, (7) substantially reduces the amplitude of the high-frequency spatial 
Fourier components of any flow variable 7. Consequently, 7 can be more accurately termed 
the large-scale part of 7. At this point, it should be mentioned that as a result of the 
defining properties of G, it follows that 


d7__ d7 dl _ dJ 

dt dt ’ dxk dxk ’ 


( 10 ) 


The turbulent fields are decomposed as follows, based on Favre filtering, 

7=7+7' ( 11 ) 

where the Favre filter 

7= y (12) 

is defined in an analogous manner to the Favre time average which has been of use in the 
more traditional studies of compressible turbulent flows (Hinze 1975). However, contrary 
to the more traditional Favre time averaging, 


7^7 

(13) 

in general, and hence 


7' ^ 0. 

(14) 

The direct filtering of the continuity equation (l) yields 


dp d(pv k ) _ 
dt dx k 

(15) 


t A Gaussian filter is adopted for the calculations in this study 
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where we have used (10) and (12). Likewise, a direct filtering of the momentum equation 
yields 



d{pv k ) d(pv k vi) _ dp do kl dr kl 

dt dxi dx k dxi dx t 

(16) 

where 

p = pRT 

(17) 



(18) 

and 

Tkl = -p(VkVl - V k Vl + VfcVj + v|vjk + v£v|) 

(19) 

is the subgrid-scale stress tensor. The subgrid-scale stress tensor 

can be decomposed as 

follows, 

T =L + C + R 

(20) 

where 

L k i = ~p{y k vi - v k vi) 

(21) 


C kl = -p{v’ k vi + v\v k ) 

(22) 


R kl = -pv k v\ 

(23) 


are respectively, the subgrid-scale Leonard, cross, and Reynolds stresses based on Favre 
filtering. From (21), it is clear that the Leonard stress can be calculated directly and does 
not need to be modeled. The cross stress is modeled with the Galilean invariant scale 
similarity model 

Cki = ~p{vkVi - v k i>i ) (24) 

which is analogous to its incompressible counterpart which has been reasonably success- 
ful in the large-eddy simulation of incompressible turbulent flows (Bardina, Ferziger and 
Reynolds 1983, Speziale 1985). The subgrid-scale Reynolds stress tensor is separated into 
deviatoric and isotropic parts, respectively, as follows: 

R = £)R + /R (25) 

where 

D R kl = -p(vM - ^vjvjtu), (26) 

and 

iRki = -^pv'MSkt. (27) 

Here, the deviatoric part of the subgrid-scale Reynolds stress tensor, j?R , is modeled using 
the compressible generalization of the Smagorinsky model that is given by 

dRh = 2C R pA 2 1 1 1 ^ 2 (S k i — -S mm 6 kl ) (28) 
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where 


c- _ 1 l dik . d ^\ 
kl 2 dxi dxk 

II 8 = SfnnSmn 


(29) 

(30) 


(i.e., S is the Favre filtered rate of strain tensor while 7/j is its second invariant) and Cr is 
the compressible Smagorinsky constant. The isotropic part of the subgrid-scale Reynolds 
stress tensor, /R , is modeled by 


iRki = —^CipA 2 IIg6ki 


(31) 


where Cj is a dimensionless constant. Equation (31) can, for the most part, be obtained by 
making a turbulence production equals dissipation hypothesis (Yoshizawa 1986). Hence, 
the overall subgrid-scale stress model we propose takes the form 


tu = ~p(vkVi - v k vi) + 2 C R pA 2 II l - ,2 (S kt - | S mm S k i ) 

- 2 -Cj-pA 2 II- s 6 M . (32) 

In the incompressible limit, the deviatoric part of (32) (which is all that is of consequence in 
an incompressible flow since the isotropic part can be absorbed into the pressure) reduces 
to the linear combination model 


— D T ki/ P = r>(vjkt > x - Vkvij — 2 Cr A 2 IIj 2 Ski (33) 

of Bardina, Ferziger, and Reynolds (1983) where the Bardina constant is one in order 
to satisfy Galilean invariance (Speziale 1985). This follows from (32) since, for constant 
values of p, 

v = V, S mm = S mm = 0 (34) 


A direct filtering of the energy equation yields the filtered form 
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c i a(itf) + 
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is the subgrid-scale heat flux. The subgrid-scale heat flux can be decomposed similarly to 
the decomposition of the subgrid-scale stresses. This leads to 


Q = Q« + 


(37) 
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where 


QW 

= C p p{v k T - vT) 

(38) 


= c p />(tjr + Cf ') 

(39) 

qW 

= C p pv[T> 

(40) 


are the Leonard, cross, and Reynolds heat fluxes. Analogous to the modeling of the cross 
stresses, the cross heat flux is modeled using the scale similarity format 

Qk ] = C p p{v k T - ~v k f) (41) 


which is Galilean invariant. The Reynolds heat flux is modeled with the usual gradient 
transport format as follows (cf. Eidson 1985) 

<5* = -C T A J /jf J2I (42) 

where Ct is a dimensionless constant such that 

Ct = ~C K (43) 

given that Pr T is the turbulent Prandtl number. Of course, the Leonard heat flux can be 
calculated directly. Hence, the overall model for the subgrid-scale heat flux we propose is 

Qk = C p p{(hT - v~ k T ) - C T A 2 1 if |£], (44) 

which is obtained by combining (38), (41) and (42). 


At this point, some comments need to be made concerning the viscous terms on the 
right-hand side of (15) and the pressure gradient-velocity and viscous dissipation terms 
which appear on the right-hand side of (35). For high Reynolds number turbulent flows 
that are sufficiently far from solid boundaries and are not in the hypersonic regime (the 
only cases that are considered for this initial study), the viscous terms described above can 
be neglected. The pressure gradient-velocity correlation can be written in the alternative 
form 


v k 
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is not yet closed. This pressure dilatation term is extremely difficult to model and not much 
success has been achieved in dealing with it. Fortunately, there exist many turbulent flows 
where the precise details of its structure, as documented in the turbulence literature, are 
not important. Consequently, at this stage in the development of compressible large-eddy 
simulations, we propose to defilter T and V -v and to calculate (46) directly. Since the large 
scales (which contain most of the turbulent energy) are likely to be the main contributor to 
the pressure dilatation term, this technique should be relatively satisfactory despite the fact 
that high frequency information is lost by defiltering. In an analogous manner, the viscous 
terms on the right-hand side of (16) and (35) can be obtained by defiltering in turbulent 
flows where the Reynolds numbers are not extremely high or where the turbulence is not 
in close proximity to solid boundaries. 

The turbulence model proposed herein is complete once values for the constants Cr, 
Ci and Ct are obtained. This will be accomplished using the results of direct numerical 
simulations of compressible isotropic turbulence. 


3. Numerical Method 


Our direct simulations of compressible turbulence are based on Eqs. (l)-(3), with the 
time derivative in the energy equation written solely in terms of the pressure. In order 
to alleviate the severe stability limit imposed at very low Mach numbers by the acoustic 
waves, a splitting method is adopted. The first step integrates the equations 
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The constants c 0 and are the current root mean square (rms) value of the sound 
speed and the free-stream Mach number, while 7 = C p /C v , where C v is the specific heat 
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at constant volume. The free-stream Mach number is These equations are non- 

dimensionalized in terms of a length scale ( lo ), a velocity scale (uo), and a pressure scale 
(p 0 ), a reference viscosity po and a reference conductivity /Cq. The Reynolds number is given 
by Re = poUolo/po and the Prandtl number is Pr = C p po/ K o • Except for the substitution 
of -j^ for p, the dissipation function $ is given by Eq. (5). Initially, the density p 0 is 
uniform and equal to one. Both the viscosity p and the conductivity k are constant in 
this initial study. The computational domain is a cube, normalized to [0, 27 t] 3 . Periodic 
boundary conditions are imposed in all three directions. 

The spatial derivatives in these equations are approximated by a Fourier collocation 
method (see, for example, Hussaini and Zang (1987)). In each coordinate direction, N grid 
points are used: x* y = 2nj/N, for j = 0, 1, . . . , N-l. The derivative of a function 
J(x) with respect to x* is approximated by the analytic derivative of the trigonometric 
interpolant of J(x) in the direction x*. Most simulations of incompressible, homogeneous 
turbulence have used a Fourier Galerkin method. The compressible equations, however, 
contain cubic rather than quadratic nonlinearities and true Galerkin methods are relatively 
more expensive (compared with collocation methods) than they are for incompressible flow. 
The essential difference between collocation and Galerkin methods is that the former are 
subject to both truncation and aliasing errors, whereas the latter have only truncation 
errors. As discussed extensively by Canuto, et al (1987), the aliasing terms are not signif- 
icant for a well-resolved flow. However, care is needed to pose a collocation method in a 
form which ensures numerical stability. For this reason, the second term in (48) is actu- 

di v x ) } d%) d f ^ 1 

ally used in the equivalent form M — * - — b pvj — — + v* — - — — 1. As noted by Feiereisen, 

OXi OXi OX\ 

Reynolds and Ferziger (1981), when this form is employed together with a symmetric dif- 
ferencing method in space (for example Fourier collocation) , then in addition to mass, and 
momentum, energy is conserved for the ideal compressible equations (zero viscosity and 
conductivity) in the absence of time differencing (and splitting) errors. 

The second fractional step of the splitting, given by (50) - (52) contains most of the 
effects of the acoustic waves. This splitting is employed at each stage of a third-order 
Runge-Kutta method. In the simulations reported here, the second fractional step is 
integrated analytically. In Fourier space (50) -(52) become 


dp A 

— + tkimi = 0 

at 

(53) 

d rni ... « n 

at +, *' p=0 

(54) 

dp 2 

— + ic 0 kimi = 0 

(55) 


where m; = pvi and Fourier transformed quantities, which depend upon the wavenumber 
k, are denoted by a circumflex. The exact solution of these equations may be written 

p( 2 ) = pW -f L\A cos(cokA.t t ) + B sin(c 0 kAt s ) — A] (56) 

c o 
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( 57 ) 


m[ 2 ^ = mj 1 ^ — — r[A sm(c 0 fcAt,) — B cos(cokAt a ) + 5] 

CqK 

p( 2 ) = A cos(c 0 kAt a ) + B sin{c 0 At a ) (58) 

where k = |k| is the magnitude of the Fourier wavenumber, A = pM, B = i^kim^\ 
The superscript 1 denotes the result of the first fractional step of the splitting and the 
superscript 2 the results of the second fractional step. The effective time-step of the 
Runge-Kutta stage is denoted by At,. The advantage of this splitting is that the principal 
terms responsible for the acoustic waves have been isolated. Since they are treated semi- 
implicitly, one expects the time-step limitation to depend upon v+ |c— c 0 | rather than v + c. 
(Although there is also a viscous stability limit for the first fractional step, it is well below 
the advection limit in the cases of interest.) This is clearly a substantial advantage at 
low Mach numbers. Since the second fractional step is integrated analytically, it does not 
contribute to any time step limitations. If one is truly interested in all the details arising 
from the sound waves, then the time-step must be small enough to resolve the temporal 
evolution of these waves. But, if only the larger-scale sound waves are of interest, then 
this splitting method is useful. 

During the acoustic fractional step an isotropic truncation is performed: for each vari- 
able — p, pv and p, all Fourier coefficients for which 

kik t > (JV/2) 2 (59) 

are set to zero. This reduces the numerical anisotropy produced by a cubic truncation. 
Moreover, it reduces the aliasing interactions in the collocation method (Canuto et al 1987, 
Chapters 3 and 7). 

The compressible code can also be executed in a purely explicit mode. In this case no 
splitting is performed; Eqns. (47)-(52) are simply combined in the appropriate manner 
and integrated directly. 

The expected stability limit of this three-dimensional Fourier collocation method for 
the compressible Navier-Stokes equations has the form 

(60) 

where for the semi-implicit version u is the maximum value, for any velocity component, 
of jv* | + \c — c 0 | and, for the explicit version, u is the corresponding maximum of |v*| + 
|c|. For the third-order Runge-Kutta method employed here the coefficient a should be 
approximately 1. In practice we have found the explicit code to be stable provided (60) 
is satisfied with a = 1.5 with isotropic truncation. The corresponding limit for the semi- 
implicit code is roughly 0.4. 

A number of simulations have also been conducted of strictly incompressible flow. These 
were performed with a separate code which also used a Fourier collocation method, but 
for the simpler, incompressible Navier-Stokes equations. 


Ax 

At < a 
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4. Comparison with Incompressible Results 


The initial conditions for the numerical simulations were designed to mimic the exper- 
imental data of Comte Bellot and Corrsin (1971), hereafter referred to as CBC. These 
experiments were also the basis of direct simulations used by Clark, Ferziger and Reynolds 
(1979), Bardina, Ferziger and Reynolds (1983), and McMillan and Ferziger (1979) in their 
analyses of LES models for incompressible turbulence. Initial conditions are chosen to 
match CBC measurements at a non-dimensional time of 240 (cf. table 4 in Comte Bellot 
and Corrsin 1971). The computational domain is a cube of side 20cm. These parameters 
are associated with measurements taken behind a grid with a mesh spacing of one inch, 
and a mean fluid velocity of 393.7 inches/sec. The initial time in the direct simulation 
corresponds to t = 0.00254 sec in the CBC experiment. The reference length (/ o), velocity 
(Uo) and pressure (po) are respectively 20cm/27r, lcm/ sec and 1 gr/cm sec 2 . After gen- 
erating a random, divergence-free velocity profile, the kinetic energy (in Fourier space) is 
scaled to match the measured CBC energy spectrum (see Appendix). Finally, the velocities 
are scaled so that the initial rms velocity agrees with the measured values. This adjust- 
ment is typically less than 1%, which provides one measure of the uncertainty in the fit 
to the experimental data. With the chosen non-dimensionalization, the Reynolds number 
Re = UqIq/uq is 22.74 based on a reference kinematic viscosity u 0 = 0.14 cm? / sec. In ta- 
ble 1, the parameters measured by CBC at t=240 are summarized. The Taylor microscale 



CBC 

64 s- 

96 s 

128 s 


6.75 

6.75 

6.75 

6.75 

ifr(Sk) 

- 

0.0 

0.0 

0.0 

E 

~ 

68.3 

68.3 

68.3 

e 

462 

375 

432 

447 

An 

0.26 

0.28 

0.27 

0.25 

a 12 

- 

0.20 

0.19 

0.27 

Ais 

- 

0.20 

0.19 

0.26 

Rx 

38.1 

43.6 

41.3 

37.8 


Table 1: Initial conditions based on CBC experiment and Clark et al. 
(1979) calculation. Mach number is zero. 


length is defined by 


1 1/2 


= 


< vj > 

< <fe* > 


(61) 
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and the dissipation £ by 


£ = 2 fij SijSijd 3 x, 


where is the rate of strain tensor 


*-s<£ + f$- 


(62) 


(63) 


In Eq. (61), <> denote a root mean square. Its exact definition is given in the Appendix. 
The Taylor microscale Reynolds number is 


R x = 




v 


(64) 


The velocity derivative skewness and flatness tensors Sk and FI are the third and fourth 
moments of the velocity gradient and are defined by 


Sk u = 


Fl u = 


,dv k .„ 

< (§^) J > 3 ' ! 
< <^) 4 > 


(65) 


(66) 


In tables 1 and 2, only the average diagonal components of the skewness and flatness tensors 
are shown. The remaining columns list the parameters obtained from the initial conditions 
of the numerical simulations on 64 s , 96* and 128* grids. There is a 20% discrepancy 
between the dissipation obtained by CBC and the dissipation computed on the coarsest 
grid which suggests that a 64* grid has marginal resolution, at best. A 12% difference 
between the value of R\ calculated on the 64 s grid and that obtained by CBC confirms 
the need for grids finer than 64*. On a 96* grid, both the dissipation and R\ are in much 
closer agreement with CBC. Discrepancies between our results and CBC for e and R\ are 
respectively 6.5% and 7.5% on a 96* grid. On the finest grids that the direct simulations 
were performed on, the computed values of e and R\ respectively have relative errors of 
3.5% and less than 1% when compared to CBC. 

The numerical simulations were run from time t = 240 until t = 375 (in CBC units) , 
which corresponds to a non-dimensional interval of 0.1145 (in our units). Table 2 furnishes 
a comparison of the experimentally measured parameters with those from the numerical 
simulation at the final time. On the coarsest grid, the total dissipation rate that was 
calculated is still slightly below the value measured by CBC. A 96 s grid generates values 
of e consistent with CBC. 

Skewness measures the departure from symmetry. At t=0.1145, the diagonal compo- 
nents of Sk are -0.5 which agree well with the numerical results of Kerr (1985). Kerr 
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CBC 

64 3 


TOM 

<7u 

5.03 

5.18 

5.19 


i(r(Sk) 

- 

-0.42 

-0.51 


E 

38.6 

40.3 

40.4 


e 

154.4 

151.3 

154.6 

156.8 

An 

0.34 

0.33 

0.34 

0.33 

A i2 

- 

0.24 

0.24 

0.23 

A13 

- 

0.24 

0.24 

0.23 

Rx 

36.6 

37.8 

40.4 

37.2 


Table 2: Final conditions (t=0.1145) based on CBC experiment and 
Clark et al. (1979) calculation. Mach number is zero. 


studied isotropic, turbulent flow, but prevented the decay of energy by using an exterior 
energy source at the large length scales. 

As noted in the previous section, we have chosen not to de-alias the advection terms. 
In reaching this decision we drew upon the extensive evidence that has accumulated on 
aliasing effects in the last dozen years (Canuto, et al 1987., Chapters 3, 4 and 7) as well 
as some tests conducted with the incompressible isotropic turbulence code. In this code, 
de-aliasing is accomplished by applying the 2/3-rule (Canuto, et al 1987, Chapters 3 and 
7) in an isotropic fashion; e.g., the de-aliased results for a 64 3 grid are obtained by running 
the incompressible code on a 96 3 grid and applying the truncation given by (59) with 32 
in place of N/2 on the right hand side. The results are summarized in Fig. 1. Here we 
present, for 64 s , 96 s , and 128 3 grids, the energy spectra E(k) (defined in the Appendix) 
at t = 0.0586 for both aliased and de-aliased calculations. Some adverse effects of aliasing 
are apparent on the 64 3 grid, but they are only in the tail of the spectra, and they are 
already insignificant on a 96 3 grid. For the reasons outlined here, a 96 s grid was chosen as 
the standard discretization for the incompressible as well as the compressible simulations. 


5. Compressible Turbulence Results 

5.1. Direct Simulations 

For the compressible case, the initial conditions imposed on the velocity distribution re- 
main unchanged. Reynolds and microscale Reynolds numbers retain the values used in 
the incompressible simulations. The initial pressure distribution over the entire field is 
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specified. The fluctuating pressure, p/ is determined from the velocity distribution by en- 
forcing a zero initial time derivative for V • v. A Poisson equation for p/ is obtained from 
the divergence of the momentum equation after setting the time variation of V • v to zero 
(Feiereisen et al. 1981). The mean pressure, p TO , is then determined so that a prescribed 
initial mean average Mach number, M 0 , defined to be the ratio of rms fluid velocity and 
rms of the speed of sound, is achieved. An analytic expression for p m is given by 

Mq f v 2 d 3 x -I- / ^-d 3 x 
P™ yf p~ 1 d 3 x 

The initial average Mach number is specified at the outset of the direct simulations (DS) 
as an initial condition. Density is initially set to unity, while the temperature, if required, 
is derived from the equation of state. 

The Mach 0.6 case contains localized regions of supersonic flow as evidenced by tables 
3-5 and by the three-dimensional contour of Mach 1 furnished in figure 2. Nonetheless, 
the statistical properties of the flow remain largely unaffected by compressibility effects if 
the average Mach number is in the low subsonic range. This is shown in figures 3-6 which 
track the time histories of a small number of statistical variables obtained from 96 3 DS. 
The time histories for skewness (fig. 3), A (fig. 4), and total kinetic energy (fig. 5) at 
Mach numbers 0.0, 0.1, 0.4 and 0.6 are virtually superimposed on each other. 

Flatness and skewness are affected the most by compressibility effects. Figure 3 indi- 
cates that the skewness which corresponds to an isotropic turbulent state monotonically 
increases with Mach number. It is -0.50 at M=0 and has increased to -0.46 at M=0.6 . 
Before the flow has reached a state of isotropic turbulence, the time evolution of skewness 
at all Mach numbers are indistinguishable from each other. The physical system has equi- 
librated after approximately one third the total computation time. While not reaching an 
equilibrium value, it is nonetheless worthwhile to point out that the flatness parameter 
decreases by 2% as the Mach number is raised from 0.0 to 0.6 as seen in figure 6. 

Figure 5 illustrates the decay of turbulent kinetic energy (1/2 as a function of time. 
This decay is a natural consequence of viscous damping. After a brief initial increase, R\ 
continuously decreases in time, (as illustrated in figure 7), with no sign of stabilizing. On 
the other hand, A which is representative of the medium-scaled eddies, increases in time. 
This indicates that energy in the higher wavenumbers is being depleted by the molecular 
viscosity. 

Tables 3 to 5 summarize the results of direct simulations of compressible isotropic 
turbulence at Mach 0.1, 0.4 and 0.6 at t=0.1145 on three different grids. Incompressible 
results are included for comparison. On all the grids, the compressible data converges to 
the incompressible results as the Mach number is driven towards zero. As expected, the 
divergence of velocity no longer vanishes, and is now an increasing function of Mo. 

While the dissipation is approximately the same on the two finer grids, the consistently 
lower values on the coarsest grid confirm the previously stated conclusion that a 64 3 grid 
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Mo 

E 

6 

< V • v > 

BS2&EB 

<M> 

M max 

0.0 

40.26 

157.4 

0.00 

-0.424 

0.00 

0.00 

0.1 

40.82 

158.2 

0.17 

-0.440 

0.07 

0.21 

Ea 

41.09 

160.4 

— B 

-0.428 

0.28 

0.84 

E9 

41.32 

162.3 

■Hi 

-0.406 

0.43 

1.26 


Table 3: Summary of direct simulations on a 64 s grid with initial 
average Mach numbers of 0.0, 0.1 and 0.4 at t=0.1145. 


M 0 

E 

e 

< V • v > 

BSJS3B 

< M > 

Mmax 

E19 

40.35 

154.3 

0.00 


0.00 

0.00 


40.49 

155.5 

0.14 


0.07 

0.23 

m 

40.79 

157.0 

1.17 

■ 

0.28 

0.93 

E9 

41.04 

158.3 

2.82 


0.42 

1.39 


Table 4: Summary of direct simulations on a 96 3 grid with initial 
average Mach numbers of 0.0, 0.1, 0.4 and 0.6 at t=0.1145 


cannot resolve all the length scales. As a function of increasing Mach number, the trace 
of Sk increases, the dissipation decreases, while the total kinetic energy increases very 
slightly. 

The results in tables 3-5 are averages over several DS runs with different initial seeds. A 
given seed uniquely determines the initial velocity distribution, and therefore the pressure 
and temperature fields. Variations of the seed are only necessary to eliminate the statistical 
uncertainty due to the random velocity distribution. The distribution of velocity on two 
different grid sizes are different even when the initial seed is the same. 

Skewness is even more sensitive to the grid refinement than the dissipation as witnessed 
by its decrease from a value of -0.505 to one of -0.521 on 96 s and 128 3 grids respectively. 
This sensitivity might be the result of the derivative of the fluctuating velocity field which 
appears in the definition of the skewness tensor Sk. 
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M 0 

E 

6 

< V- v > 

< 5 <r(Sk) > 

<M> 

Mmaz 

0.0 


153.98 

0.00 

jmsEgm 

msM 

0.00 

0.1 


154.84 

0.14 

KB 

ISO 

0.21 

0.4 

40.78 

156.54 

1.11 

■SO 

B9 

0.86 


Table 5: Summary of direct simulations on a 128 3 grid with initial 
average Mach numbers of 0.0, 0.1 and 0.6 at t=0.1145 


5.2. Data Analysis 

Using the data generated from the previously discussed DS of compressible homogeneous 
turbulence at low Mach numbers, the proposed subgrid-scale (SGS) model is validated. 
If a good model is to be successful, it must require that the subgrid stresses calculated 
from the DS correlate well with the computed model stresses. Models express the subgrid 
scale stresses, not available to a large eddy simulation code, as a function of the large 
eddy velocities which are known. These velocities are simply the Favre-filtered velocities 
introduced earlier. The Favre-filtered velocities are calculated by filtering the resolved DS 
velocity field with a Gaussian spatial filter of width A/Ax/, where A x/ is the grid spacing 
on the fine grid. Perturbed velocity fluctuations on the fine grid are the difference between 
the fully resolved velocity and the filtered ones, 

v' = v - v. (68) 

From v' and v subgrid stresses based on DS, considered exact, are calculated. These 
include the Leonard (21), cross (22) and Reynolds (23) subgrid stresses. However, these 
subgrid-scale stresses themselves do not directly affect the evolution of the system. The 
momentum equation is only influenced by the divergence of the subgrid stresses (vector 
level). Similarly, v • V • r (scalar level) is a better representation of the dissipation terms 
in the energy equation than are the stresses. Consequently, correlations are performed on 
the tensor, vector and scalar levels. High correlations must be achieved on all levels of 
► comparisons before the SGS model can be considered useful. 

The data analysis proceeds in multiple stages. First, the exact stresses calculated from 
the DS are injected down to the coarse grid, along with with the filtered velocities. The 
modeled subgrid stresses are then calculated (excluding the model constants) on the coarse 
grid. Some variables must be filtered a second time (e.g. cross stress terms). Rather than 
calculate them on the fine grid (which is not available to the large-eddy simulation codes) , a 
Gaussian filter is applied to v on the coarse grid with a filter width of A c Ax c . Consistency 
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between the coarse and fine filter widths is achieved by insuring that 


A L = N l 
A c N c 


( 69 ) 


where Nf and N c are respectively the number of nodes along one direction on the fine and 
coarse grids. This guarantees that the filtering on the coarse and fine grid is performed over 
the same region in physical space. Henceforth, the coarse and fine filter widths are refered 
to in their non dimensional form A c and A/. Derivative evaluations on both the coarse 
and the fine grid are based on Fourier collocation. Calculations by McMillan and Ferziger 


(1979) indicate that the model constants are sensitive to the accuracy of the derivative 
evaluations. A general trend that has been observed is that the Smagorinsky constant 
is lowered when derivative quantities are evaluated more accurately. Our constants are 
therefore expected to lie in the lower range of the values obtained by McMillan (1980). 


Next, the two model constants, Cr andC/ , are calculated. Unfortunately, the constants 
can be calculated by a wide variety of algorithms, each with its own merits. Moreover, for 
each algorithm, the constants can be evaluated from tensor, vector or scalar information. 
Therefore, criteria must be established to identify the best method. A key test is that 
L + R should be Galilean invariant. To make use of this fact, an additional constant, 
Cc , is introduced as an extra factor in the subgrid cross stress model. A self consistent 
method of calculating the constants will have to produce a Cc of 1 to satisfy the Galilean 
invariance property stated above (Speziale 1985). Additional tests are performed on coarse 
grids with varying degrees of refinement which further decrease the number of choices. A 
thorough discussion of model constants is the subject of the next subsection. 


Once a single or a multiple set of model constants have been determined, the model 
subgrid stresses are calculated and correlated with the exact subgrid-scale stresses calcu- 
lated from the DS after injection onto the coarse grid. The correlations are performed for 
each type of subgrid-scale stress individually, and for the total stress (L + C + dR + /R ). 
Strong differences in the correlation coefficients relating total stresses are noticed depend- 
ing on whether or not the Leonard stresses are included. Finally, the correlations obtained 
from the proposed model are compared with the linear combination model, which has been 
shown to be the best model available for incompressible isotropic turbulence. Correlation 
coefficients are calculated based on the two pairs of constants that are obtained from the 
above considerations. The set that is finally retained corresponds to the highest levels of 
correlation of the total stress on the vector and scalar levels. These matters are treated 
more completely in a later subsection. 


To avoid a possible confusion of terminology when referring to variables being compared 
against each other, superscripts () m and () e are sometimes used. They respectively refer 
to modeled and exact (based on DS) variables (tensor, vector or scalar level). 
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5.2.1. Model Constants 


The proposed model (32) has two undetermined coefficients. The first constant, C R , is 
associated with the modeled deviatoric subgrid Reynolds stress, £>R , while the second 
constant, Cj , allows the modeled isotropic stress, /R to be adjusted. 

Although the cross stress model has no constant associated with it, it is nonetheless 
multiplied by a constant Cc . This is done in the hope of reducing the number of schemes by 
which the constants can be calculated. A good model should reproduce a cross constant of 
one to guarantee the Galilean invariance of C m . Once the constants have been determined, 
Cc is set to one and forgotten. Because the flow is isotropic, constants are expected to 
be the same for the three diagonal stress components, the three off-diagonal components 
and the three vector components. Therefore, the values presented in the tables below are 
averaged over the appropriate components. In the tables, D refers to averaged diagonal 
components, OD to average off-diagonal components, V to averaged vector components 
and S to scalar components. Similar averagings are performed for correlation coefficients. 

The three constant models (Cr , Cj and Cc ) are calculated using two techniques, each 
applied on the tensor, vector and scalar levels. One method enforces equality of the rms 
levels of the exact and the modeled stress. This is done for each individual subgrid- 
scale stress (deviatoric and isotropic parts of the subgrid-scale Reynolds stress and the 
cross stress). Hereafter this approach is referred to as RMS. The second method utilizes 
multiple linear least square regression (LSQ) between the exact (19) and the modeled (32) 
total subgrid-scale stress to determine the constants. Table 6 summarizes the results 
obtained from incompressible data. The three cases presented are identical except for 
the initial random number seed. The constants are independent of the detailed velocity 
statistics, except for Cj iS<? which varies randomly. These results are based on a vector 
level comparison between the modeled and exact stresses. Both RMS and LSQ produce 
C c close to unity as required. Unfortunately this prevents a rational choice from being 
made between the two approaches. A more complete set of LSQ constants is presented in 
table 7. They are computed at Mach numbers of 0.0, 0.1, 0.4 and 0.6 on a tensor, vector and 
scalar level. Computations are performed on a coarse grid of 16 3 . Confirming the results 
of table 6, Cj has strong variations as a function of Mach number on the vector level when 
the random seed is changed. Additionally, Cj is the subject of random fluctuations at the 
scalar level. On the tensor level, Cj is quite independent of Mach number and the random 
seed. A possible interpretation of the variation is that the isotropic stresses do not have 
much influence on the total stress correlations. That this is indeed the case is explained in 
the section on correlations. The remaining constants are stable. That is to say, they are 
very little influenced by either Mach number or initial conditions. 

At first glance, Cc is near unity at both the scalar and vector levels. However whereas 
on the vector level, the constant remains within 4% of unity for all Mach numbers, this is 
not the case on the scalar level where Cc is a decreasing function of Mach number. This 
trend is present in the data generated from both random seeds. Although not presented 


17 



LSQ 

RMS 

Seed 


r Cj 


C R 

Cj 

Cc 

1 

Esta 

0.0070 

isa 

0.023 

0.036 

1.03 

2 

m 

0.0024 

EH 

0.022 

0.034 

1.04 

3 

0.012 

0.0055 

1.03 

0.022 

0.033 

1.03 


Table 6: Model constants calculated by LSQ and RMS between ex- 
act and modeled total stresses (L + C + R ) . Results are based on 
three identical incompressible simulations except for the random ini- 
tial velocity distributions. Calculations are on the vector level on a 
16 3 grid. 


here, the RMS cross stress model constant is also near unity when calculated based on 
vector level stresses. Therefore, a preferred method for the determination of the model 
constants is still not possible. 


5.2.2. Filter Width and Grid Coarseness 

Confirmation of the numerical evidence presented by McMillan and Ferziger (1979) that 
A c = 2 is the best filter width is given in table 8 (Mo = 0.1). The criterion used to 
determine the validity of the filter width is that Cc RMS must remain close to unity on the 
vector level. Only when A e = 2 is Cc near one. Similar results hold for LSQ constants. 
The constants also vary with respect to the coarse grid on which the LES is to be per- 
formed. Table 9 summarizes the model LSQ and RMS constants evaluated from modeled 
stresses calculated on 16 s and 32 3 grids on the vector level. The data (Mo = 0.1) shows 
that Cr varies by 30% when the grid size on which the modeled stresses are calculated 
ranges from 16 3 to 32 3 . On the other hand, large eddy simulations using finite difference 
algorithms might be performed on grids as large as 128 s . Unless a subgrid scale model is 
found which produces constants independent of the coarse grid size, LES simulations will 
run the risk of producing the wrong results. Perhaps a more complicated dependence of 
the modeled subgrid Reynolds stresses on A is required. 

The best model constants will produce the highest correlations between the modeled 
and exact subgrid stresses on all levels (tensor, vector and scalar). Based on the previous 
discussion, only constants calculated on a vector level are adequate because they produce 
a Cc of unity. Unfortunately, a clear cut choice between RMS and LSQ constants cannot 
be made because Cc (on the vector level) is nearly one in both cases. Rather than make 
an arbitrary choice, both sets of constants are considered when correlating exact against 
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modeled stresses. For reference, the constants used henceforth are 

C R Lsq = 0.012 Ci LSQ = 0.0066 

C r rms = 0.023 Cj RMS = 0.035. 


(70) 

(71) 



Seed 1 

Seed 2 

Mo 


0.0 

0.1 

0.4 

0.6 

0.0 

0.6 


D 

0.12 

0.12 

0.12 

0.12 

0.11 

0.11 

Ci 

V 

— 

0.0070 

0.0067 

0.0063 

0.0023 

0.0023 


S 

— 

.0034 

.00012 

-0.25 

.0017 

.0003 


D 

1.32 

1.32 

1.32 

1.31 

1.32 

1.31 

C c 

V 

1.04 

1.04 

1.02 

1.00 

1.04 

1.00 


S 

1.00 

1.00 

0.95 

0.873 

0.971 

0.934 


D 

0.018 

0.018 

0.018 

0.018 

0.016 

0.015 

C R 

V 

0.012 

0.012 

0.012 

0.012 

0.012 

0.012 


S 

0.015 

0.015 

0.015 

0.014 

0.015 

0.015 


Table 7: LSQ model constants. Filter widths are A/ = 12 and A c = 2 



LSQ 

RMS 

El 

6 

12 

mm 

6 

12 

mm 

H 

1 

2 

■■ 

1 

2 

mm 

Cr 

0.007 

0.012 





Ci 

-0.02 

0.002 

0.014 




C c 


1.03 

1.33 

0.82 

1.03 

1.13 


Table 8: LSQ and RMS model coefficients between exact and modeled 
stresses on the vector level at Af 0 = 0.1. Results are obtained with 
fine filter widths of 6, 12 and 24 while maintaining the proper ratio 
of 6 between fine and coarse widths The coarse grid is 16 3 . 


5.2.3. Correlations 


Correlation coefficients have long been a preferred diagnostic tool for estimating the re- 
liability of the modeled stresses. However, the subgrid stresses have been separated into 
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Table 9: LSQ and RMS model constants based on subgrid stresses 
evaluated on 16 s and 32 s coarse grids 


various components (L, C,R), each modeled separately and although the correlation of 
any of these components against their models might be excellent, it is still possible for the 
correlation of the exact total stress against the modeled total stress to be less impressive. 
Such is the case, for instance, when two stress components have opposite signs and par- 
tially cancel each other out. As a final note before the specific correlation coefficients are 
presented, one must always be attentive to the actual relationship between the exact and 
the modeled variable, even when the correlation coefficient is relatively high (say, above 
70%). A correlation coefficient as high as 70% may not be as good as it seems. For exam- 
ple, the correlation between the parabola y = x 2 and the constant y = 1 in the interval 
[0, l] is 75% although they are uncorrelated. As a consequence, correlations are deemed 
good if the correlation coefficient is above, at least the 90% level. 

For convenience, the compressible subgrid scale model is restated here: 


D Rii = 2C s /SA ! CSk.Su )' 12 (Sii - &W<i) 

(72) 

/&/ = -|C,?A ! (SuSu) 

(73) 

Cij = -7f{viVj - ViVj) 

(74) 


The correlation coefficients presented in table 10 between exact and modeled /R , 
dR and C, are independent of the model constants. The correlation coefficients are in- 
sensitive to the average Mach number variation, except for the isotropic stress correla- 
tion coefficients which exhibit higher values for the larger Mach numbers. Although the 
isotropic stresses enjoy a correlation above 80% on the tensor level, the correlation coeffi- I 

cients plummet below 15% on the vector and scalar levels (we will see later that these low 
correlations are not of great consequence). This indicates that the isotropic modeling is 
poor. 

In most of the literature on subgrid-scale models, the Leonard stress has been omitted 
from the total stress correlations on the grounds that it is calculated exactly ( Bardina, 

Ferziger, and Reynolds 1983, McMillan 1980, McMillan and Ferziger 1979). However, it has | 
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M 0 


0.1 

0.4 

0.6 


D 

31 

31 

31 

£)R 

OD 

26 

26 

25 


V 

22 

22 

22 


s 

45 

45 

45 


D 

89 

89 

89 

c 

OD 

91 

91 

91 


V 

80 

80 

80 


s 

75 

74 

72 


D 

84 

85 

85 

I R 

V 

13 

15 

18 


S 

12 

14 

18 


Table 10: Correlations between exact and modeled stresses. Shown 
are the correlations for £>R , jR , and C at Mach 0.1, 0.4 and 0.6. 


recently been shown by Speziale (1985) that the combination L + C is Galilean invariant, 
while the cross stress is not. Therefore we feel that correlations of the total stress should 
include the Leonard stress. Table 11 summarizes the correlation coefficients between the 
total stress including and excluding the Leonard stress. Results are presented at Mach 0 
and Mach 0.6. When the Leonard stresses are left out, correlation coefficients similar to 
Bardina’s are obtained on all levels, even though Bardina did not include the contribution 
of the isotropic stresses. However, the inclusion of L decreases the correlations at the 
vector level by approximately 30%. Table 11 substantiates that the correlation coefficients 
(with and without L) are nearly independent of the initial average Mach number. 

Correlation coefficients between (C + R) e and various combinations of the modeled 
stresses using constants calculated by LSQ and RMS are summarized in table 12. The 
second column indicates the model against which the total stress is being compared. Al- 
though at first glance RMS based constants perform better at a tensor level when the 
total modeled stress is considered, at the vector and scalar levels, the trend is reversed. 
Correlations at the vector and scalar level are higher by 4% using the LSQ constants. The 
only cases where this trend is absent is when either R or /R is removed from the modeled 
total stress in which case the choice of constants has a negligible effect on the correlation 
coefficients. 

When the constants are selected based on LSQ, table 12 indicates that the correlations 
on all levels are highest when all modeled components are included. The inclusion of 
/R only affects correlations on the tensor level. However, the dynamic evolution of the 
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M 0 = 

0.0 

Mo = 

0.6 


L+C+R 

C+R 

L+C+R 

C+R 

D 

93 

82 

93 

81 

OD 

80 

85 

79 

84 

V 

46 

72 

46 

71 

S 

56 

73 

56 

74 


Table 11: Comparison of correlation coefficients of the exact total 
stress versus its model. The modeled terms are computed on a 16 3 
coarse grid. Each case is presented with and without the inclusion 
of the Leonard subgrid-scale stress terms (calculated exactly). Both 
the incompressible and the M=0.6 case are shown to illustrate the 
weak influence of Mach number on the correlation coefficients. Model 
coefficients are Cr = 0.012 and Cj = 0.0066. 


large scale velocities only brings into play the stresses on the vector and scalar level. 
Therefore the coefficients which produce the maximum correlations of total stress on these 
two levels should be chosen. This leads to an optimum choice of 

C R = 0.012 

Cj = 0.0066 (75) 

calculated by least squares fit of the total stress on the vector level. 

Conversion of Cr to the standard form currently used in incompressible LES * produces 
a Smagorinsky constant of 0.092. McMillan (1979) obtained a value of Cs — 0.10 when 
spectral collocation derivative computations were employed. This constant corresponds to 
a scalar level evaluation based on RMS. On the vector level, McMillan calculated a higher 
Smagorinsky constant of 0.13. This value can be obtained from the present data by using 
£>R rms instead of D R LS<i . However as shown above, the correlations of the total subgrid 
stress would be lower. 

Yoshizawa (1986) recently performed a direct interaction approximation calculation 
based on a model for slightly compressible turbulent flows. His model for the isotropic 
component of the subgrid Reynolds stress is identical to ours. However he calculates 
Ci = 0.177, to be compared to our value of 0.0066 s . Yoshizawa also obtains a Smagorinsky 

*In a number of reports (McMillan and Ferziger 1979, Bardina et al. 1983), the subgrid Reynolds stress 
model is proportional to Cs 2 . In these cases, the relationship between Cs and Cr is Cr = y/2 Cs 2 . 

8 This comparison must be viewed cautiously since Yoshizawa (1986) did not base his model on Favre- 
filtered fields 
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Least squares 

RMS 

C R 

=0.0122, Cj 

=0.0066 

Cr 

=0.023 , C/ 

=0.035 

Exact 

Model 

D 

OD 

7 

S 

D 

OD 

V 

S 

C+R 

C 

78 

82 

68 

62 

78 

82 

68 

62 


C+/R 

81 

82 

69 

62 

89 

82 

67 

61 


C+x>R 

78 

84 

71 

70 

77 

83 

69 

69 


C+R 

81 

84 

71 

70 

88 

83 

67 

67 


Table 12: Correlation of the exact total stress (C+R ) with various 
models. The modeled stresses are defined in equations (72)-(74) 



constant Cs = 0.17 which is consistent both with incompressible results and the results 
in this paper. 

Results presented in table 12 indicated that the isotropic stresses do not noticeably 
influence the modeled total stress on the vector and scalar level. When the Navier-Stokes 
equations are written in the incompressible form, the isotropic stress and the trace of 
the cross stress can be included with the pressure to form a new total pressure variable. 
Although this cannot be done for compressible flows because the pressure satisfies an 
evolution equation of its own, it seems plausible that the isotropic stress should exert 
very little influence if its rms value is much smaller than the rms of pressure and if the 
same relation holds true between the gradients of these quantities. At Mach 0.6, the rms 
value of the pressure and pressure gradient are 283 and 92, respectively, as opposed to 
the corresponding values of 7.1 and 4.4 for the isotropic Reynolds subgrid-scale stress. 
These values verify the fact that this term does not play an important role. At lower 
Mach numbers, the ratio between the rms isotropic stresses and the rms pressure is even 
greater, and tends to infinity in the incompressible limit. At M=0 it is well known that 
the isotropic stresses have no influence on the velocity distributions. 

Total stress correlations obtained with the compressible model are compared with re- 
sults calculated with the linear combination model and are summarized in table 13. Data 
is based on the LSQ model constants. The most meaningful comparison is at M = 0.0 
where the exact, and modeled cross and Reynolds stress are isotropized. The first two lines 
of the table show a very good agreement between Bardina’s results and ours. However, 
the third line, which includes jC in the modeled and exact total subgrid stress, generates 
a correlation coefficient of 70% on the scalar level, or 10% higher than that obtained by 
Bardina. 

Initial tests of the subgrid-scale heat flux model produced a turbulent Prandtl number 


23 



Exact 

Model 

Stress 

Vector 

Scalar 

Bardina 

£>C + £>R 

£>C +z?R 

83 

74 

60 

M 0 = 0.0 

L >C + £)R 

£>C +z)R 

86 

77 

60 

M 0 = 0.6 

c+ d r 

c+ d r 

85 

73 

71 


Table 13: Comparison of the linear combination model (row 1) with 
two different combinations of subgrid-scale stresses at Mach 0.0 and 
0.6 (rows 2 and 3). Correlations were obtained on a coarse grid of 
16 s . 


of 


Pvt — 0.4 


(76) 


which is a little low compared to the expected value of approximately 0.7. This may be 
due to the fact that the temperature gradients are not very large in this problem, and 
hence, these numerical results may not be well suited for the accurate calculation of Pr T . 
This issue will be addressed in a future study. 
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6. Conclusion 


New subgrid-scale models for compressible turbulence have been developed and tested 
against the results of direct numerical simulations of compressible isotropic turbulence. 
These compressible subgrid-scale models, which were based on the Favre-filtered equations 
of motion for an ideal gas, contain three dimensionless constants and reduce to the linear 
combination model in the incompressible limit. The subgrid-scale stress model constants 
were found to assume the values of Cr — 0.012 and Cj = 0.0066 and those gave rise to 
correlations between the exact and modeled stresses that were above 70% on the tensor, 
vector and scalar levels - a correlation which compares favorably with those obtained in the 
previous subgrid-scale modeling of incompressible turbulent flows. Another encouraging 
feature lies in the fact that these constants and their associated correlations were found 
to be comparatively insensitive to the Mach number in the range 0 < Mo < 0.6. While 
the modeling of the isotropic part of the subgrid-scale Reynolds stress is relatively poor, it 
appears that it does not play a significant role in the dynamics of the flow for the range of 
Mach numbers considered in this study. Likewise, the turbulent Prandtl number obtained 
herein (Pry = 0.4) appears to be somewhat inaccurate. We believe that this problem can 

; alleviated by conducting correlations on the basis of turbulence computations which 

✓olve more intense temperature gradients (an effort which must await future research) . 

Future large-eddy simulations of compressible turbulent flows are planned based on 
the subgrid -scale models developed in this study. Undoubtedly, as such studies progress, 
additional refinements in these models will be introduced. Nevertheless, we believe that 
the essential foundations for the large-eddy simulation of compressible turbulent flows have 
been established in this study. This could have a profound impact on the future analysis 
of supersonic and hypersonic flows in aerodynamics. 
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Appendix 

Initial energy spectrum 

Both the incompressible and compressible direct simulations impose a specified energy 
spectrum on the initial random velocity distribution. For comparative purposes, the data 
was obtained from tabular data found in Comte-Bellot and Corrsin (1971). They tabulate 
the function Fn(fc) which is related to the energy spectrum E(k) by 

(Al) 

Unfortunately the data is noisy, so a least squares fit is performed on log Eu, expressed as 
a fourth order polynomial in log k. The final form obtained for E n is 

log (Eu) = 2.64359 - 0.72602(log £) -0.32585 (log A:) 2 

+0.03525(log k ) 3 - 0.02344(log fc) 4 . (A2) 


Calculation of model constants 


Before correlating the total exact subgrid stress against its model, the model constants 
must be determined. There are many ways of accomplishing this among which two are 
retained. The total modeled stress is written as a linear combination of modeled terms 

C<r?: 

T m = f2 C ir, m (A3) 

«=i 

while the exact total stress is simply 

r e = X>/. (A4) 

1=1 


The unknown constants to be determined are the C,. The first approach adopted is to 
calculate the root mean square of the pairs C,r t m and r/ and equate them for each value 
of i. The constants thus take the values 


(J T rn 



This method is referred to by RMS in the text. 
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A least square method applied to the total stress as a whole is an alternative approach. 
In this case, the norm 

ll>- m -r'||’ = ||E(r?-C i r')|| ! (A6) 

t=l 

is minimized with respect to the coefficients C,. This gives rise to a linear system in the 
coefficients which can be solved by direct methods if the number of constants is not too 
large. For the subgrid model considered in the text, n = 3. 

Definitions 

For reference purposes, several statistical definitions are provided here. All variables are 
defined on a three-dimension grid and are subscripted by a single index * for convenience. 
The average of a function is 


<?>= (at) 

As a function of the average, the rms of J is 

a 7 = y/< (/- < 7 >) 2 >. (A8) 

Correlation coefficients are fundamental in evaluating subgrid-scale models. The correla- 
tion coefficient between two functions Jand £is 

C(7,$) = <TS> . (A9) 

OjOg 
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Fig. 5 Time history of one third the trace of kinetic energy for M 0 of 0.0, 0.1, 0.4 and 0.6. 
Direct simulation was performed on a 96 3 grid. 
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Fig. 6 Time history of one third the trace of FI for M 0 of 0.0, 0.1, 0.4 and 0.6. Direct 
simulation was performed on a 96 s grid. 
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Fig. 8 Scatterplot of the x component of modeled (C + R ) m versus (C + R ) c on the vector 
level. Results are based on a direct simulation on a 96 3 grid and a data analysis on 
a 16 3 grid. 
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